home *** CD-ROM | disk | FTP | other *** search
/ Practical Algorithms for Image Analysis / Practical Algorithms for Image Analysis.iso / CH_4.4 / XFM / XFM.C < prev    next >
C/C++ Source or Header  |  1999-09-11  |  6KB  |  282 lines
  1. /* 
  2.  * xfm.c
  3.  * 
  4.  * Practical Algorithms for Image Analysis
  5.  * 
  6.  * Copyright (c) 1997, 1998, 1999 MLMSoftwareGroup, LLC
  7.  */
  8.  
  9. /*
  10.  * XF(ast)M(oments)
  11.  *
  12.  * calculate moments, M, of a gray scale image up to m33,
  13.  * based on a filter array, Y, extracted from the imput image,
  14.  * and a mask CM:
  15.  *
  16.  *           M = CM*Y*CM'; CM': transposed CM, *: matrix mult.
  17.  *
  18.  * ref:      M. Hatamian, 
  19.  * IEEE Trans. on Acoust., Speech and Signal Proc., 
  20.  *    ASSP-34, 546 - 553, 1986;
  21.  * Proc of Digital Signal Proc. - 87,
  22.  *    V. Cappellini & A. G. Constantinides, Eds.
  23.  *    Elsevier Sci. Publ. B. V. (North Holland), 1987
  24.  *
  25.  * version based on M. Hatamian's code implementing recursion 
  26.  * to obtain Y-matrix;
  27.  *
  28.  */
  29.  
  30. #include <stdio.h>
  31. #include <stdlib.h>
  32. #include <math.h>
  33. #include "xfm.h"
  34.  
  35. #define    ON        1
  36. #define    OFF        0
  37.  
  38. #define    RESET        ON               /* reset graphics screen */
  39. #define    NO_RESET    OFF
  40.  
  41.  
  42. #undef    DEBUG
  43.  
  44.  
  45. extern char *optarg;
  46. extern int optind, opterr;
  47. extern short tiffInput;         /* flag=0 if no ImageIn to set tags; else =1 */
  48.  
  49.  
  50. int FILTER_TEST = OFF;
  51. int WRITE_FILE = OFF;
  52.  
  53.  
  54.  
  55. void
  56. fail_alloc (str, code)
  57.      char *str;
  58.      int code;
  59. {
  60.   printf ("\n...memory alloc for %s failed\n", str);
  61.   exit (code);
  62. }
  63.  
  64.  
  65. /*
  66.  * usage of routine
  67.  */
  68. void
  69. usage (char *progname)
  70. {
  71.   progname = last_bs (progname);
  72.   printf ("USAGE: %s inimg [-d] [-w file] [-t] [-L]\n", progname);
  73.   printf ("\n%s iteratively evaluates moments of a convex shape up to 3rd order,\n\n", progname);
  74.   printf ("ARGUMENTS:\n");
  75.   printf ("       inimg: input image filename (TIF)\n\n");
  76.   printf ("OPTIONS:\n");
  77.   printf ("          -d: default mode: evaluate moments of binary region\n");
  78.   printf ("          -w: write file (.mdt) to disk\n");
  79.   printf ("          -t: generate test filter array\n");
  80.   printf ("          -L: print Software License for this module\n");
  81.  
  82.   exit (1);
  83. }
  84.  
  85. void
  86. main (argc, argv)
  87.      int argc;
  88.      char *argv[];
  89. {
  90.   int i, j;
  91.   int nr, nc;
  92.  
  93.   int jmin, imin, jmax, imax;
  94.   int left_x, right_x;
  95.  
  96.   float m00;
  97.   float mu00, mu11, mu02, mu20;
  98.   float xc, yc, rg, dent;
  99.  
  100.   int xdim, ydim;
  101.   float *y[MDIM], *m[MDIM];
  102.  
  103.   char *buf;
  104.   int i_arg;
  105.   FILE *file;
  106.   Image *imgIn;                 /* input image */
  107.  
  108.  
  109. /* 
  110.  * cmd line options
  111.  */
  112.   static char *optstring = "dtw:L";
  113.  
  114.  
  115. /*
  116.  * parse command line
  117.  */
  118.   optind = 2;
  119.   opterr = ON;                  /* give error messages */
  120.  
  121.   if (argc < 2)
  122.     usage (argv[0]);
  123.  
  124.   while ((i_arg = getopt (argc, argv, optstring)) != EOF) {
  125.     switch (i_arg) {
  126.     case 'd':
  127.       printf ("...option %c: default mode\n", i_arg);
  128.       FILTER_TEST = OFF;
  129.       break;
  130.     case 't':
  131.       printf ("...option %c: generate test filter array\n", i_arg);
  132.       FILTER_TEST = ON;
  133.       break;
  134.     case 'w':
  135.       printf ("...option %c: write file %s to disk\n",
  136.               i_arg, buf = optarg);
  137.  
  138.       if ((file = fopen (buf, "w")) == NULL) {
  139.         puts ("\n...could not open output file");
  140.         exit (1);
  141.       }
  142.       WRITE_FILE = ON;
  143.       break;
  144.     case 'L':
  145.       print_sos_lic ();
  146.       exit (0);
  147.     default:
  148.       printf ("...unknown command line argument encountered\n");
  149.       usage (argv[0]);
  150.       exit (1);
  151.       break;
  152.     }
  153.   }
  154.  
  155. /*
  156.  * get the input image
  157.  */
  158.   imgIn = ImageIn (argv[1]);
  159.   if (imgIn->bps == 8 && imgIn->spp == 3) {
  160.     printf ("Got RGB image!!!\n");
  161.     fprintf (stderr, "Can only work with Grayscale or Binary TIFF files!!!\n");
  162.     exit (1);
  163.   }
  164.   if (imgIn->bps != 1) {
  165.     printf ("Got Grayscale image!!!\n");
  166.     fprintf (stderr, "Can only work with Binary TIFF files!!!\n");
  167.     exit (1);
  168.   }
  169.  
  170.   /* reset tiffInput so that we write a grayscale file (i.e tags are not copied) */
  171.   tiffInput = 0;
  172.  
  173.   jmin = imin = 0;
  174.   jmax = imgIn->width;
  175.   imax = imgIn->height;
  176.  
  177.   left_x = jmin;
  178.   right_x = jmax - 1;
  179.   nc = jmax - jmin;
  180.   nr = imax - imin;
  181.  
  182. #ifdef DEBUG
  183.   printf ("\n...nc = %d, nr = %d\n", nc, nr);
  184. #endif
  185.  
  186.  
  187. /*
  188.  *  zero outer border of image
  189.  */
  190.   zero_border (imgIn, 2);
  191.  
  192.  
  193. /*
  194.  * memory allocation for filter and moment arrays
  195.  */
  196.   xdim = ydim = MDIM;
  197.   for (j = 0; j < ydim; j++) {
  198.     if ((y[j] = (float *) calloc (xdim, sizeof (float))) == NULL)
  199.         fail_alloc ("y", 1);
  200.   }
  201.   for (j = 0; j < ydim; j++) {
  202.     if ((m[j] = (float *) calloc (xdim, sizeof (float))) == NULL)
  203.         fail_alloc ("m", 1);
  204.   }
  205.  
  206. /*
  207.  * begin moment evaluation
  208.  */
  209.   if (FILTER_TEST == ON) {      /* test filter array */
  210.     printf ("...initialize a test filter array\n\n");
  211.     for (j = 0; j < ydim; j++) {
  212.       for (i = 0; i < xdim; i++) {
  213.         if ((j == 0) || (i == 0))
  214.           *(y[j] + i) = (float) 1.0;
  215.         if ((j == 1) && (i == 1))
  216.           *(y[j] + i) = (float) 1.0;
  217.         if ((j == 1) && (i == 2))
  218.           *(y[j] + i) = (float) 1.0;
  219.         if ((j == 2) && (i == 1))
  220.           *(y[j] + i) = (float) 1.0;
  221.         if ((j == 2) && (i == 2))
  222.           *(y[j] + i) = (float) 1.0;
  223.  
  224.         printf (" %2.0f ", *(y[j] + i));
  225.         if ((i + 1) % 4 == 0)
  226.           printf ("\n");
  227.       }
  228.       printf ("\n");
  229.     }
  230.   }
  231.   else                          /*  compute filter output from given image */
  232.     filter_recursion (left_x, right_x, nr, nc, y, xdim, ydim, imgIn);
  233.  
  234.  
  235. /*
  236.  * evaluate moments by matrix multiplication
  237.  */
  238.   printf ("\n...matrix multiplication to obtain moments\n");
  239.  
  240.   eval_mom (y, m, MDIM);
  241.  
  242.  
  243. /*
  244.  *  compute centroid and central moments
  245.  */
  246.   central_moments (m, &mu00, &mu11, &mu02, &mu20, &xc, &yc, &rg, &dent);
  247.  
  248.  
  249. /*
  250.  *  print area-normalized and central moments
  251.  */
  252.   printf ("\n...moments:\n");
  253.   m00 = *(m[0] + 0);
  254.   for (i = 0; i < (int) MDIM; i++) {
  255.     for (j = 0; j < (int) MDIM; j++) {
  256.       *(m[i] + j) /= m00;
  257.       printf ("  m[%d][%d]/m00 = %f\n", i, j, *(m[i] + j));
  258.     }
  259.   }
  260.   printf ("  %f\n", m00);
  261.   printf ("\n...central moments:\n");
  262.   printf ("...mu00 = %f,    mu11  = %f\n", mu00, mu11);
  263.   printf ("...mu02 = %f,     mu20 = %f\n", mu02, mu20);
  264.   printf ("...rg = %f\n", rg);
  265.   printf ("...dent = %f\n", dent);
  266.  
  267.  
  268. /*
  269.  * write data file of format ( .mdt)
  270.  */
  271.   if (WRITE_FILE == ON) {
  272.     write_mdt_file (file, m, (int) MDIM, m00, mu00,
  273.                     mu11, mu02, mu20, xc, yc, rg, dent);
  274.     fclose (file);
  275.   }
  276.  
  277.   for (j = 0; j < ydim; j++) {
  278.     free (y[j]);
  279.     free (m[j]);
  280.   }
  281. }
  282.